The Spread Optimization Problem 
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We analyze the irreversible dynamical process corresponding to the linear-threshold model of 
influence spread over a network, and propose an efficient algorithm to solve the inverse problem, 
namely that of finding an optimal initial condition that generates to a desired final state of the 
dynamics. This is a challenging problem which is crucial in many contexts, from systemic risk 
analysis to the design of viral marketing campaigns. Our approach is based on the cavity-method of 
statistical physics. We compare our algorithm to standard techniques (based on Monte Carlo and 
Linear Programming) and show that it has an excellent performance. 



INTRODUCTION 

Much effort has been devoted over the past forty years 
to investigate the properties of irreversible dynamical sys- 
tems, such as their relaxation or the dependence of the 
equilibrium state on the initial condition. Recently the 
interest extended also outside physics, with the applica- 
tion of statistical physics methods to problems arising in 
economics, biology and computer science. One limitation 
of the approaches taken so far is that they provide ways 
to explain and predict the observed phenomena, but no 
instrument to control them. This aspect may be unusual 
in physics, but it is central in the applied sciences. In 
engineering for instance, the main goal is that of design- 
ing the structure of a system or the initial conditions of a 
dynamical process in such a way to favor/prevent the oc- 
currence of some phenomenon. It follows that one of the 
challenges in the interdisciplinary research on complex 
systems should be that of developing efficient methods 
and algorithms to control and optimize their dynamical 
behavior. This is particularly important when these sys- 
tems are large networked structures with many degrees 
of freedom and complex interaction patterns [T] . 

The simplest form of dynamics taking place on net- 
works is probably that of irreversible spreading processes, 
in which nodes "activate" depending on the state of their 
neighbors. Despite the straightforward mathematical for- 
mulation, spreading processes find application in a va- 
riety of relevant practical problems. For instance, one 
could be interested in activating the whole graph at time 
T using a minimum number of "seeds" . This is a typi- 
cal problem faced in the design of viral marketing cam- 
paigns, whose goal is that of targeting those individuals 
that guarantee the maximum spread of the advertisement 
throughout a social network [H [3] . A similar optimiza- 
tion process can be used to identify which nodes are most 
sensitive for the propagation of failures in a power grid or 
liquidity shocks in interbank networks, providing a tool 
to assess and prevent systemic risk [3]. 

In this Letter, we show that the cavity method and 



message-passing algorithms can be applied to efficiently 
solve spread optimization problems on large graphs, pro- 
viding a framework to study also more complex inverse 
dynamical problems. Our results demonstrate that the 
performance of the optimization depends on the nature 
of the collective dynamics and, in particular, on the pres- 
ence of cooperative effects. Moreover, the dynamical tra- 
jectories selected under optimization are characterized by 
peculiar statistical properties, such as a very slow pace 
with large fluctuations of the activation times. We con- 
clude showing the relation between our Belief Propaga- 
tion (BP) based solution of the inverse dynamics and 
recursive equations recently used to study the dynamics 
of other irreversible processes on networks. 



THE SPREAD OPTIMIZATION PROBLEM 

We consider a linear-threshold model (LTM) !.2j in 
which the nodes of a graph G = {V,E} can determin- 
istically switch from a quiescent state {xi = 0) to an 
active one {xi — 1) during a discrete-time dynamics. 
Each directed edge (j,j) G E \s endowed with a weight 
Wij G R+, while the nodes are given a set of thresholds 
{9i E R+, i e V}. A quiescent node i activates when the 
total weighted influence of its active neighbors exceeds 
its threshold, i.e. it follows the discrete-time parallel dy- 
namics 



r-t + i 



1 if xl 

otherwise 



1 or EiGSz "'j^s;* > 



(1) 



We shall call "seeds" the nodes active at t — 0. 

The direct dynamics for this model has been recently 
solved in the context of the zero-temperature random 
field Ising model [7] . Here we adopt a different approach 
allowing us to study an associated optimization problem 
that focuses on large deviations from the typical case 
analyzed in these works. We shall highlight later the 
relations between our approach and the non-optimized 
case. 



Because of irreversibility, a trajectory {x°, . . . , x'^} can 
be fully parametrized by a configuration of activation 
times t = {ii, . . . , tj^}, with U G {0,1,2,..., T, 00} rep- 
resenting the time at which Xi becomes equal to 1, and 
conventionally defining i^ = 00 iff cc* = for any t <T, 
an arbitrarily defined stopping time. Given an assign- 
ment of the seeds S = {i : ti — 0} , the dynamics is fully 
determined for i ^ S" by the set of relations ti = 4>i{{tj}) 
with j e di and 4({ij}) = minEjeai«)jafe<t]>0. * 
(we conventionally assume the minimum of an empty 
set to be 00). Defining '^i{U,{tj)) = l[ti = Q] + 
1 [ti — (/)i({tj})], the dynamical constraints can be simply 
expressed by the equations ^j(ij,{ij}) — 1 for every i. 
To optimize over the assignment of the seeds, we weight 
each dynamical trajectory with an energetic term (?(t) = 
Z^i ^i{ii) where <Si{ti) = ^qI [U = 0] + el [U = 00] and 
Ci e]0, 1] is the cost to be selected as a seed. The trade-off 
between the total cost of the seeds and the final number 
of active nodes is controlled by the chemical potentials 
/J,, e S R. The resulting optimization problem is com- 
putationally hard even to approximate in the worst case 



THE BP SOLUTION 



In order to eliminate systematic short loops due to the 
fact that the constraints ^i,\E'j share the variables ti^tj 
for (i,j) G E, we employ a dual factor graph in which 
variable nodes are associated to edges (i, j) € E and rep- 
resent the pairs of times [ti , tj ) while the factor nodes 
contain three terms: 1) the dynamical constraints, 2) 
the constraint imposing that ti is the same for all edges 
incident on z, and 3) the contribution from i to the en- 
ergy (see Supplemental Material, henceforth SM). Us- 
ing this representation, the optimization problem can be 
analyzed on locally tree-like graphs by means of the cav- 
ity method and solved using message-passing algorithms. 
In fact, the (cavity) marginal probability Hii{ti,ti) that 
nodes i and £ activate at times ti and t^ in absence of ^^ 
satisfies the following BP equations 

H,,(i„i,)cxe-''^-(*')^*,(t„{t,})ni/,,(tj,t,) (2) 

with j £ di\ i. Given a solution of (|2|, the marginal 
probability of {ti,tj) is Pij{t^,tj) ex Hij(ti,tj)Hj^{tj,t^). 
Eqs. (I2]) allow to study the statistics of the dynam- 
ical trajectories (e.g. entropies or distribution of ac- 
tivation times). For /3 — >■ cxj, with a proper rescal- 
ing of the messages, they give the Max- Sum (MS) al- 
gorithm, which can be used to find explicit solutions 
to the spread maximization problem (see SM). Some 
kind of approximation is often needed in practice to ob- 
tain efficient updates. For small connectivity, continuous 
weights and thresholds must be discretized, in order to 





FIG. 1. Example of low energy configurations found by MS 
equations with reinforcement in 2D L x L lattices with open 
boundaries. Colors correspond to activations times (from 
time zero for seed in red to colder colors). Left: hexagonal 
lattice (k = 6) with 6 = 5,T = W, L = 21. Every hexagon 
is adjacent to at most one hexagon of a colder color. Right: 
rectangular lattice {k = 4) with 6 = 3,T = 10, L = 20 



have only a finite spectrum of possible values for the sum 

J2jedi'^JiM^j < *]• ^^ '-^^^ ® *^^ integer value of the 
discretized threshold. We notice that the BP-messages 
do not depend on all T^ combinations of times {ti^t^), 
but only on (ij,sign(i,j — if + 1)). It follows that the 
time complexity per node and per iteration of BP and 
MS is 0{TQ'^\di\\di-l\). By exploiting the associativity 
of the convolution needed to perform the local update 
rule (see SM), the overall time complexity per iteration 
can be reduced to scale as 0{TQ'^\E\). In the case of 
BP, the convolution can then be performed with Fast 
Fourier Transform, reducing further the complexity to 
0{T<d'^\E\). While it is clear that increasing the reso- 
lution on the weights and thresholds could be computa- 
tionally very expensive, the computational complexity of 
message-passing is still very low compared to other opti- 
mization algorithms such as exhaustive search and Monte 
Carlo methods. On the other hand, for large connectiv- 
ity (e.g. fully connected graphs), continuous weights can 
be kept, and the update can be simplified by using the 
Central Limit Theorem. 

Fig. [1] shows some examples of low energy config- 
urations obtained by MS on two-dimensional lattices. 
Note that the cavity approximation is not expected to 
be asymptotically exact on lattices, but that the config- 
urations found have reasonably low energy. 



COMPARISON WITH OTHER METHODS 

We tested on random graphs two alternative strategies 
to solve the optimization problem of minimizing the en- 
ergy S": a Linear/Integer Programming (L/IP) approach 
and Simulated Annealing (SA). We tried two different 
L/IP formulations: one closely related to our represen- 
tation based on x\ and the same dynamical constraints 
'^i and a second one based on [TS] , solving both with the 
CPLEX software. In both cases performances were ex- 
tremely poor (becoming unfeasible for graphs over a few 
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FIG. 2. Comparison between the performances of SA and 
MS algorithms on random regular graphs of degree K — 5, 
thresholds 6 = 4 and different sizes N. We computed the 
quantity Apo ~ Po /po ~ 1 (averaged over 10 realiza- 
tion of the graph) that represents the rescaled difference be- 
tween the minimum density of seeds po required to activate 
the whole graph obtained using SA with exponential cool- 
ing scheme from /? = 0.5to/3 = 10'^ and the corresponding 
values Po^ computed using the MS algorithm. The results 
for different system sizes A'^ are plot as function of the total 
number of global updates M employed in the SA. Points are 
slightly shifted horizontally from the correct value at A^ = 100 
for clarity. The running time of SA for each instance with 
N = 1000,M = 102 400 was around 73 hours. For compar- 
ison, the same instances were solved by MS in less than 1 
minute on the same machine. 



dozen nodes) and will not be reported in detail. With 
SA vi^c find that the time to reach the MS solution scales 
poorly with the size (see Figure [2]) when the effective con- 
figurational space is sufficiently large. A detailed investi- 
gation of the performances of the algorithms depending 
on the choice of parameters (costs, weights, etc.) is be- 
yond the scope of this Letter. 



COLLECTIVE DYNAMICS UNDER 
OPTIMIZATION 

For random initial conditions (i.e. i.i.d. seeds proba- 
bility), there exists a density of seeds Po<l above which 
the whole graph activates during the dynamics. Depend- 
ing on the weights and the thresholds, the complete ac- 
tivation can occur continuously in po or abruptly, as a 
result of some macroscopic cooperative mechanism. In 
the following, we show that the effects of the optimiza- 
tion are very different in case of continuous or discon- 
tinuous activation processes. For the sake of simplicity 
we consider instances of regular random graphs of degree 
K, with uniform weights Wij = 1, y{i,j), uniform thresh- 
olds 6i = 9, \/i with 9 G {1, . . . ,K—1}, and uniform costs 
Ci = 1, Vi (an especially difficult case for BP/MS due to 



the high degeneracy) . The dynamics of the LTM is now 
reminiscent of bootstrap or fc-core percolation processes 

[6] . For random initial conditions, the collective be- 
havior of the system depends on the values of 9. For 

1 < 9 < K — 1, the whole graph activates (pr = 1) 
abruptly at a finite density pg, whereas for 9 = 1 and 
9 = K — I there is a continuous percolation-like transi- 
tion. 

We studied the ensemble of solutions of the problem in 
the single link approximation using population dynamics. 
In the completely homogeneous case (with Ci = l,Vi), 
the population simplifies to a single message. In general, 
MS and BP results prove that the optimization is effec- 
tive when the activation process occurs continuously at 
e = 0. In particular for e > a finite interval of values 
of < p < Pc(e) exists in which full-spread solutions are 
achievable at po < Po- On the other hand, when the acti- 
vation is discontinuous at e = 0, the single-link equations 
do not converge in the region of /i < Pc(e) where one can 
expect an improvement in the value of po required for full 
activation. 

FigJ3]A displays the curves pT vs. po ioi 9 = 2 and 
K = 3 with and without optimization. The results ob- 
tained in the absence of optimization (e = 0) correspond 
to the solution of the dynamics starting from random 
initial conditions in which x° = 1 with probability pq. 
While the activation is continuous for e = 0, increasing 
e > we observe the emergence of a gap with the upper 
branch that optimizes px- The gap is just a manifesta- 
tion of a first-order transition that is evident in Fig[3j3, 
where we plot the curves po{p.) for e = (black line) and 
0.4 (red dashed line). Following the two branches of so- 
lutions across the transition point we find a coexistence 
region in which the upper branch becomes metastable 
(see blue lines in FigjSB) and hysteresis phenomena can 
be observed. The convergence of the BP equations slows 
down in the optimized region where convergence time 
grows very fast by increasing p along the upper branch. 
On the contrary, for decreasing p in the lower branch the 
convergence is fast until the transition point where the 
equations do not converge anymore. The convergence 
problems are amplified by both increasing the inverse 
temperature /3 and/or the duration time T and can be re- 
duced using a damped update rule or (large) populations 
of messages. 

In Fig|4] we report the activation time statistics of the 
dynamical process. For random initial conditions, the ac- 
tivation probability P{t) of a node decays exponentially 
with the time t as it usually occurs in spreading pro- 
cesses. Interestingly, in the region in which optimization 
is effective, the spreading process proceeds at a slower 
pace. For e > and increasing p from negative to posi- 
tive values, the dynamics slows down and P{t) develops 
a power-law behavior that persists until p ~ pc- 

The results seem different in the case 9 = 2 and isT = 4, 
in which the complete activation occurs abruptly even 
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FIG. 3. A) Parametric plot pr v.s.po obtained solving Eqs. H 
in the single-link approximation on regular random graphs of 
degree K = 3, for threshold 6 = 2, duration T = 20 and 
e = 0,0.1,0.4,1. The vertical arrow indicates the minimum 
density of seeds (po ~ 0.253) necessary for the total activa- 
tion obtained by MS on finite graphs of size \V\ = 10, 000. B) 
Curves po(m) for e = (black full line) and 0.4 (red dashed 
line). The blue curves are obtained following the upper and 
lower branches of solution across the transition. C) Activa- 
tion time probability P{t) obtained computing the total BP 
marginals in a dynamics of duration T = 100, for e = 0.4 and 
difi'erent values of /i. 



for e = 0. As reported in Fig|4]'\., the region of complete 
activation does not widen when increasing e from zero, 
while the onset of discontinuity moves towards smaller 
values of po and the BP equations do not converge in a 
finite interval of /i values (see Fig|4l3). The analysis of 
the activation statistics for e = reveals important in- 
formation on the dynamics close to the abrupt transition. 
In Fig|4t! we plot P{t) for /i increasing in the optimized 
region: the non-monotonicity of P{t) for fj, > means 
that after a period of time in which the dynamics is slow 
and correlations are built throughout the graph, there is 
a sudden activation of a finite fraction of nodes that is 
reflected into the appearance of a bump in P{t) at large 
times (that diverge for /i -^ fic)- 

The above results show that, in case of continuous ac- 
tivation processes, some trajectories can be selected in 
order to optimize the spread of the process. But how 
rare are the optimal trajectories? The answer is given 
in FigJ5] (top) where we plot the entropy s of the ini- 
tial conditions that lead to a full spread. This quan- 
tity is compared with the entropy associated to a ran- 
dom choice of initial conditions with fixed density of ac- 
tive nodes po (binomial sampling without optimization). 
When the curves deviate from the binomial the probabil- 
ity of choosing randomly an optimal set of seeds becomes 
exponentially small (inset in Figjsl). In the infinite sys- 
tem, this event is governed by a zero-one law. The results 
for N = 50, 100 are obtained using a generalization of the 
cavity method that allows to fix a global constraint (the 
number of seeds) by introducing an additional set of mes- 



FIG. 4. Node activation probability P{t) obtained computing 
the total BP marginals in the single-link approximation on 
regular random graphs of degree K = 3 (top) and K = 4 
(bottom), threshold 6 = 2 with T = 100, e = 0.5 and different 
values of /i. 



sages that flow over a spanning tree superimposed on the 
original graph (see 13" Appendix B), the one for A^ = 30 
by explicit enumeration. 



RELATIONS TO EXISTING STUDIES 

It is now clear that irreversible processes on graphs 
can be optimized using BP messages that depend on 
pairs of activation times for neighboring nodes. However, 
recent works concerning the zero-temperature dynamics 
of the random field Ising model [7j and the susceptible- 
infected model of epidemic spreading [S] suggest that in 
the absence of optimization, when the seeds are randomly 
drawn, the dynamics can be correctly analyzed using only 
single-time cavity marginals. In fact, for e = 0, the two- 
times formalism can be reduced to a set of cavity equa- 
tions for Xie{t) = P'T-iixl = 1}, i.e. the probability that 
node i is active at time t in the absence of the neighbor L 
To show this, it is convenient to introduce a new represen- 
tation of the dynamic rule. For each directed edge («,^), 
the variable tu represents the time at which node i would 
activate in the absence of node (., and for i ^ S", it satisfies 
the iterative equation Ui = niinj^^_^g.^^„^..i[t^..<t]>e. t. As 
it happened for the activation times {ii}, also the equa- 
tions for the cavity-times admit a unique solution for a 
given choice of S, which is in one-to-one correspondence 
with the solution of the single ti model. In order to opti- 
mize the trajectories and average over the initial condi- 
tions, we introduce the messages Huitu^tu) that satisfy 
a set of BP equations (see SM). When e = 0, the hypoth- 
esis that the messages {Hi({tii,tu)} do not depend on 
the backward cavity times {tn} is self-consistently sat- 
isfied, and the BP equations can be easily reduced to 
single-time quantities Hnitu). Integrating the resulting 
equations for tn < t and averaging the contribution of 
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FIG. 5. Entropy s of the solutions of the full spread problem 
on regular random graphs of degree K = 3 and threshold 
9 = 2, for N = 30, 50, 100, oo vs. seed density po- For each 
A^ the upper line corresponds to the binomial distribution (i.e. 
entropy of seeds in the absence of optimization) and the lower 
one to the entropy of fully spreading seeds. Inset: Probability 
P of randomly selecting a fully spreading set of seeds for the 
same set of parameters, i.e. the ratio of the corresponding 
two curves in the main plot. 



all nodes, it is straightforward to obtain an evolution 
equation for the density p{t) of active nodes at time t. 
Although the details of the model are different, the final 
single-time equations are equivalent to those derived in 
Ref. [T (see also Ref.[5^). We remark that the two-times 
joint probability in Eqj2]is more than a technical artifact; 
on the contrary, it is quite crucial to allow information 
to flow backwards in time when optimizing over the final 
state (e > 0). 

It is worth noting that our formulation in terms of 
activation times t is a direct simplification, seemingly 
valid only for microscopically irreversible processes, of 
the general dynamic cavity formulation recently proposed 
to study non-equilibrium dynamics on sparse graphs flOl- 
[T^ (see SM). It could be in principle possible to relax 
the assumption of complete irreversibility, using proba- 
bility distributions over time-dependent paths, but the 
optimization of fully reversible dynamics is currently nu- 
merically unfeasible and it requires new separate ideas to 
overcome current computational limitations. 

Conclusions. - The study of inverse dynamical prob- 
lems on large graphs is a new frontier for the applica- 
tion of message-passing algorithms and the approach pre- 
sented here is directly applicable to a series of real- world 
problems. Nevertheless, our methods could be improved 
in two main directions. First, by including the opti- 
mization over some stochastic parameter [9], providing 
a formulation that is closer to a class of stochastic opti- 
mization problems that received a lot of attention in the 
computer science community 2J. 
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SUPPLEMENTAL MATERIAL 



Derivation of the Max-Sum equations 



The Max-Sum (MS) messages hii{ti, ti) are defined in terms of tire Belief Propagation messages Huiti^ tg) as 

hit{ti, ti) = lim - log Hii{ti,ti) + Cu 

/3— >-oo p 



where the additive constant Cu is such that 



miixhaiU,te) = 0. 



Taking the /3 — > oo limit of the BP equation (2) in the article we obtain 



max 

{tj,jedi\i} 



E %fe'0) 



jedi\ 



fJ-Ci + Cyii 



hie{ti,ti) = < 



max 

{tj,jedi\e} s.t: 
Efcea, «'fc>l[*fc<t.-l]<e> 






(-'ii 



if U = 0, 



if < U < T, 



max 

{tj,jedi\i} s.t: 



jedi\e 



1 + Cig ii ti = 00. 



(3) 
(4) 

(5a) 
(5b) 

(5c) 



These equations are derived in a straightforward way, but their computational complexity is exponential in the 
connectivity of the node being considered. We shall now see how to reduce them to a form which can be computed 
efficiently. The same method can be easily extended to the BP equations. 



Simplification of the update equations 



The maximization in (5a I is unconstrained, and it reduces to 



jedz\e 



max hji {tj , 0) 



HCi + Ci- 



(6) 



which is trivial to compute. 

The constrained maximisations in ( |5b[ [5c| can be computed efficiently by iteration. To compute ( 5b ) we introduce 
(with an obvious simplification of notation) : 



Qi,...,r{Oi,02,t) = max 

{tl,...,tr} s.t: 

^j = l....,r"'jMtj<t-'i-]=dl, 

Ej = i,...V"'jl[*i<*-ll=«2 



j = l, ■■■,'■ 



max • 

t^ 



; {hritr, t) + Ql....,r-l {dl - Wrl [tr < t - 1] , 02 

with the first element directly obtained from the definition: 

Qi{0i,d2,t) ^ max hi{ti,t) 



Url[tr <t-l],t)} 



max 

il s.t: 

u;il[ti<t-l]=ei 

iuil[ti<t-l]=e2 



( max hi(ti,t) if 6*1 = 6I2 = 0, 

ti>t-i 



(7) 
(8) 
(9) 



= < 



hi{t - 1, t) if 6*1 = wi and 02 = 0, 

max hiitiA) ii 01 = 02 ^ Wi, 
ti<t-i 

—00 otherwise. 



(10) 



The convolution of two Q's is given by 

Ql,...,r{dl,02,t) = 



max 



U-^ ,U-^ ^Uvy ,[/r, S.L. 



{Ql,...,s(^l,^2:^) + Qs+l,...,r{()li02,t)} . 

Once Q{9i,62,t) is computed from the convolution of all the incoming edges, we can easily compute 



M {61,62,1)= max { V Hj{tj,t) 

Ej = i ,u.jl[t3<t-l]>ei, V3-^:-,r 

Ej=i,...',,-«'ji[*j<*-i]<e2 
max Q{9[,6'2,t) 

{9[,e'2}^s.t: 
6^^6i,92<C92 



in terms of which ( 5b ) gives: 



htiiU, ti) = M{6, - wul [ti <U-l],6i- Wi^l [ti < i, - 1] , U) + C,j 
To compute (|5c|) we introduce 



^l.^.r^ 



max 

{tl,...,tr} s.t: 
Ej = i ,t"il[tj<T] = 



^ hj{tj,co) 



max |/lr(ir, C^) + 'S'l....,r-l(^ ~ Wrl [tr < ^] ) } 



with first clement 



'niax{/ii(r, cx)), /ii(oo,oo)} if = 
Si(9) = I i^^^ hi{ti,oo) a 6 = wi 

00 otherwise 



and generic convolution 



Once S is computed for all the incoming edges, we can compute 



max ASu...A()') + Ss+i,...An}- 

S.t: t7 -j-u^ ^ ^u 



J2 ^j'^^o 



in terms of which (5c I becomes: 



P(6) = max X / ,1,1,1,1,1 

{il,. ..,*.} s.t: I .J-^ ■' ■' 

i:,=i,...,,,«'jife<T]<e U— i'---'- 

= max S(6') 
8' s.t: e'<e 



h^eioo, te) = P{6^ - wul [h < T]) - 1 + Cu- 



if < t, < T . 



(11) 



(12) 
(13) 

(14) 

(15) 
(16) 



(17) 



(18) 



(19) 
(20) 

(21) 



Notice that P{6) = M{0,6,t) with t = 00, so we can compute it in the same manner as M{6i,62,t) provided we 
extend the range of values for t so as to include i = 00. 



Simplification of the messages 

The update equations can be further simplified by noticing that in (Is]) the dependence of ha on ti is almost trivial: 
for fixed ti, hu is independent on ti if ti = 0, it only depends on sign (7^ — {ti — 1)) € {—1, 0, 1} if < i^ < T , and it 




FIG. 6. Dual factor graph representation for the Spread problem 



only depends on 1 [ti < T] if i^ = oo. We then introduce: 

hii{ti,a) - 



hiiiti, ti — 2) if ii > 1 and a = 

h-ie(ti, ti — 1) if ii > and cr = 1 

hie(t,,ti) ifCT = 2 

— oo otherwise 



which can be inverted as 



hie{ti,ti) — hii[ti,a{ti,ti)) 
&{ti,ti) = 1 + sign (tg - (tj - 1)) 



and in terms of which ([8|), (10), (16) and ( |17[ ) become respectively: 

2, t) = max < hr(tr, (j{tr, t)) + Ql,...,r-1 (^1 ~ Wrl[^r < ^ ^ 1] 



Q 



l,....rl"i,P2, 



2,t)^ I 



max hi(ti,a{ti,t)^ ii 9i = 62 = 

ii>i — 1 

hit -1,2) 



max hi(ti,2) 
ti<t~i 



if 9i = wi and 92—0 
if 9i — 92 = wi 

-00 otherwise 

Si,...,r{0) = max |/i^(tr, 2) + Si^,,,^r-i{0 ~ w,.l[i,. < T]) | 



S'i(6') = <( max hi{ti,2 
' ti<rf-i 



-00 



|/ii(d-l,2),i?i(d,2)| if 61 = 

if 61 = wi 
otherwise 



and ^, ([U]) and ([2l]) become 



Wr^ltr < t 



1]^0} 



hie{ti,a) 



Yl 'aiaxhj,{tj,a{tj,0)) 

j€di\i '- ' 

M(6'i - wgi, 9i - Wii, ti) + Cu 
M{9i - ■wu-,9i,ti) + di 
M{9^,9,,ti) + Ca 

P{9{) -l + Cu 

—00 



^Ci + Cii if ti = and a — 2, 

if 1 < t, < T and cr = 0, 
if < tj < T and a = I, 
if < tj < T and cr = 2, 
if <i = 00 and ct = 0, 
if ti = 00 and ct > 0, 
otherwise. 



(22) 



(23) 
(24) 



(25) 
(26) 

(27) 
(28) 



(29a) 

(29b) 
(29c) 
(29d) 
(29e) 
(29f) 
(29g) 



Time complexity of the updates 



A naive implementation of 29a|[29g would require 0{Td{d~ 1)6'*) operations for a vertex of degree d, but this can 
reduced substantially to 0{Td&^) as follows. A factor d—1 can be saved by pre-computing the convolution Qi.... i 



be reduced substantially to 0{TdO"') as follows. A factor d—1 can be saved by pre-computing the convolution Q 
and Qi^...^d for each i — l,...,d using 2d convolution operations, and then computing "cavity" Qi^....i-i.i+i,...d as a 
convolution of Qi,...,i-i and Qi+i,...,d- Moreover, a factor Q^ can be saved by employing FFT to compute convolutions 
in the case of BP. 

For the case of MS, the problem can be relaxed as follows in order to reduce complexity. We consider the relaxed 
constraint ^I^^ = l[ti = 0] + l[ii > 4>{{tj}j£i)] instead of ^i. The formal consequence is that the equivalent of Eq. [s] 
has only one inequality, and thus the complexity of the max-sum convolution becomes 0{TdQ^). Notice that in the 
case of e = 0, this is equivalent to considering the following dynamics: at each time t, select a subset of the nodes and 
update them using the standard dynamical rule. In the optimization case e > unfortunately this is no longer true. 



BP equations on Regular Random Graphs 

We consider here the ensemble of random regular graphs of degree K on which we define for simplicity a completely 
homogeneous LTM with uniform weights Wij — 1 V(i,j) G E, uniform thresholds 6i — 6 \/i G V, and uniform costs 
Ci — 1 Wi € V . Because of the homogeneity of the problem, in the infinite size limit (N — > cx)) one can consider a 
single-link approximation, obtaining the following system of nonlinear equations: 

H{t,s)o. E ^",^' , pf-^'"-"°mrg(t~l,t)"° for 0<t<T 

n_-t-n++no — /^— 1 

n_<6l-l[s<t-l] 

e-l[s<i-l]<n_+no 

H(0,s)(xe-%^-i 

H{oo,s)oce-^' Y^ (^^^^[H{T,oo)+H(oo,oo)]'^-^-''-m2^ (30) 

„_<e_i_i[s<T] \ ~ ^ 

where we defined the cumulative messages pt — X]s>t ^i^T^') and rrit — X]s<t-i H{s,t). The normalization constant 
is just the sum of all messages. The system of equations could be further simplified from 0{T'^) messages to 0{T) by 
exploiting the fact that H{t,s) = H{t, sign{t — s + 1)). As this reduction requires the introduction of slightly more 
complex normalization constant, for simplicity here we do not consider it explicitly. 

In general one should study the behavior of Eqs|30| varying //,e,/3 and T for any given assignment of K and 6. 
Here we consider only two exemplar cases: 1) K = 3,9 = 2, and 2) K = 4,6 = 2. The behavior in the (e, /i) -plane is 
studied at fixed T = 20 and ;9 = 1. 



BP equations in the cavity-times representation 

In the cavity-times representation of the dynamics, a variable ta represents the time at which a node i would 
activate in the absence of node £. Let us define: 

/» {{tk}kedi} =ramh:J2 Wk^l [tk < A > 9, i (31) 

I kedi J 

f^J {{tk}ke^^\,) - min J t : ^ Wk^l [h <t]>9A (32) 
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when i is not a seed, tn satisfies the iterative equations tu — fa I {tk}k(:di\i ) ■ ^^^ messages Ha {tu, tu) satisfy the 
following BP equations 



Hie {tie,tei) ex _ 



Assuming Hk-^^ {tki,y) = Hk->.i [tki) we get 






We will show that in the case e = 0, the equation becomes equivalent to those derived for the zero-temperature 
random-field Ising model by Ohta and Sasa f^. For that purpose, we need to consider time-cumulative quantities, 
such as the probability xui't) = X^t <t^if(iif) that node i is active at time t in the absence of node £. Let us first 
define the following sequence of increasing sets Ut 



Vt = I {tk^} ■■ J2 ^fe'^ t*''^' < t] > 61, A J2 '^k^l [tk^ < t - 1] < 9, 

[ keai\e kedi\e 

Ut = I {tk^} ■■ J2 ^fe'^ [*'»' <t]>Oi 
[ keai\i 

Ut+i = UtU Vt+i 

Ut n Vt+i = 

0<t'<t 



Yl ^[*''M{tki}kediv)) 



0<t'<t 



With this definitions, let us analyze the evolution equations Eqj33]in the e = case. It is easy to see that under 
the hypothesis that all Hki{tki, tik) on the rhs. do not depend on tki and e = 0, then the dependence on tu disapears. 
This implies that the hypothesis is self-consistent and will be verified at every iteration if it is verified at the initial 



one (e.g. if the messages have uniform initialization). Now Eq. 33 becomes 
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Let us compute the time-cumulative quantities: 

Xidt) = X! ^^^ (^^'^^ 
tii<t 



oc 



E I E s{taJu{{tk^},e^^v))+e-''\ n ^'^^ (*''-) 



E 1 



{tfcilfcgaAj 



kedi\t 



_[_[ Hki{tki 
k<£di\l 



E E n 5{xkA[tk^<t])l 



{*ki}keai\i {2:fe=0,l} fcsOA 



E Wfc,l [tfc, < t] > 



kedi\i 



E 1 

E 1 
{2:^=04} 

E 1 

{2^=0,1} 



E ^fc«^fc ^ 

kGdiy 

E ^fcja;fc > 

E "^kiXk > 
kediV 






]^ < a^fc E ■^'^^ (^'^*) + (1 ~ 2;^) E ^fe^ (^fe* 



tki>t 






Denoting by a, = e '^ /{I + e '^) the probabihty of choosing i as a seed and the initial conditions are Xui^) — ^i ^* 
and fixing the normalization factor (1 + e^^)^^ = 1 — a^, we obtain 

X^tit+l) = a, + {l-ai) E ^^ XkWk^>0.:\ W [XkXk^{t) + {I - Xk){l - Xk^{t))] (34) 

{a:fc=0,l} k£di\t kedi\l 

From which we can compute finally the density of active nodes p{t) — -^^J^i Pii^)- Apart from the obvious differences 
in the details of the dynamical update, the single-link version of Eq. ^3 (assuming all Xki identical) is substantially 
equivalent to Eq. (4) in Ref . |7] . 



Relation to the Dynamic Cavity Equations 

A general method to study non-equilibrium dynamical processes on sparse graphs has been recently introduced by 
several authors under the name of dynamic cavity method |10H12j . It is easy to show that in the absence of optimization 
(e = 0), the belief-propagation equations presented in the main text can be viewed as a simplified version, valid only 
for microscopically irreversible processes, of the dynamic cavity equations. To this end, we adopt a formulation similar 
to that used by Neri and Bolle [TTl by considering the path probability P(x*|h*) of a trajectory x* = {x°, . . . ,x*} 
with X* = {xq, . . . , x]^} in the presence of an external field h* = {h°, . . . , h*}. On a cavity graph, in which node i and 
its interactions are removed, the path probability can be written as 



Pi^%') - Pi^M%' + ny n ^ l^tl^^'^f^] Poi^°) 



(35) 



where U/^> is an auxiliary external field acting over the neighbors of i and introduced to keep track of the directed 
infiuence of i over them alo ng the dynamics. For simplicity we have taken a factorized distribution over the initial 



conditions. By summing Eq 35 over all possible trajectories of all nodes j ^ i, we get the local marginal ^i(a;*|h*-l-u^s) 
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that satisfies the equation 

t 

t 



(36) 



where y^^ is the joint trajectory of the neighbors of i up to time i, and the last passage conies from the fact that the 
dynamics is parallel and the transition probability for variable i at time t depends only on the states of neighbors at 
time t — 1. Exploiting the tree- like assumption, we can express the path probability distribution over the neighbors 
of i in factorized form 



jedi 



(37) 



In the present case, we have hf = 6i and m^j — WijX^ for all times s > 0, while the transition probability is given 
by the deterministic update rule as follows 






E '^'''^k ^ < 



.kedi 






(38) 
(39) 



For a given directed edge (i,j), thus the local field can be univocally parametrized in terms of the variable in 
the removed node, leading to the following set of recursive equations for the cavity marginals Pi^i{xl\h^ + u*^) = 

-1 i— >i!(3ii: 3i£Ji 



p,^,{xix',)^poix°) Y. n^Ni^'"';^'^] n p,^.i^v^^~'^ 

jeai\i 



(40) 



{^'r'},ea,v''=^ 



A path probability is now a vector of 0(4*+^) variables, but further reduction is possible because the dynamics is 
microscopically irreversible. The full sequence of t + 1 binary values taken by a;* = {a;^, . . . , a:*} can be encoded in a 
single integer ti — {0, ■ ■ ■ ,t} representing the time at which the variable i flips from to 1 (with the convention U = 
for a seed). With this parametrization, the dynamic cavity equations 40 take the form 



Pi^e{ti,ti) Gcpo{St^fl) Y 1 

{tj}jeoi\e 



^ Wjafe <U-l]<9^ 



jedi 



^wj-afe <u]> 



jed'i. 



-pciA't 



E 1 



{*j}jeai\f 



Y^jiMtj <U~l]< 
jedi 



jedi 



jedi\e 

n p,^.,{t„u) (41) 

jeai\e 



where we represented the initial factorized distribution in terms of weights over the seeds, i.e. po{St-fl) — e^^'^^^'i-" /{1 + 
g~tJ.ciSt-,oy n jg easy to check that Eq 41 corresponds to the e = case of the BP equations derived in the main text 
for the cavity messages Hi^i{ti,ti). From the above derivation it is also evident that optimization could be included 
by introducing in the dynamic cavity equation a local energy term £ = J2i^ii^i) — ~J2iJ2s<t^^x*.i in order to 
account for optimization over the trajectories. In the activation time representation the local energy term correctly 
becomes £i = eSt^, oo- 



